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Abstract 

We present a new Vlasov code for collisionless plasmas in the nonrelativistic regime. 
A Darwin approximation is used for suppressing electromagnetic vacuum modes. 
The spatial integration is based on an extension of the flux-conservative scheme, 
introduced by Filbet et al. [23]. Performance and accuracy is demonstrated by com- 
paring it to a standard finite differences scheme for two test cases, including a Harris 
sheet magnetic reconnection scenario. This comparison suggests that the presented 
scheme is a promising alternative to finite difference schemes. 
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1 Introduction 



For many interesting problems in plasma physics a deeper understanding can 
only be gained with the help of kinetic plasma models. In contrast to macro- 
scopic models where only the streaming velocity, pressure, temperature and 
other fluid quantities are considered, the kinetic description deals with the 
particle distribution function /(x, v). Here x is the position in space and v is 
the three dimensional velocity. The equation describing the evolution of the 
six dimensional distribution function in phase space for a collisionless system 
is the Vlasov equation. The complexity of kinetic plasma descriptions arises 
because of the electromagnetic fields appearing in the Vlasov equation. These 
depend on the charge density and current density in the plasma which are 
in turn given by the moments of the distribution function. Thus the system 
presents itself as a strongly coupled non-linear system of partial differential 
equations. Analytical solutions of this system are only possible for a very re- 
stricted number of special problems. 
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Among the numerous problems, where collisionless kinetic plasma simula- 
tions are important, we name only two topics, which are important in both 
laboratory- and astrophysical settings. These are the thin current sheet region 
in collisionless magnetic reconnection processes and the structure of collision- 
less shocks. The Harris sheet model [1] gives a kinetic equilibrium of a current 
sheet separating two regions of different magnetic field. While this can be for- 
mulated analytically, the reconnection process, which is initiated by small per- 
turbations from the Harris equilibrium, is not yet fully understood. Also, there 
has been tremendous progress in the last 10 years, both on the fluid level (see 
[2,3,4,5,6,7]) and using kinetic simulations (see [8,9,10,11,12,13,14,15,16,17] 
for Particle in Cell simulations and [18,19] for Vlasov simulations). Still many 
questions remain to be solved, such as: the spontaneous onset of reconnection, 
three dimensional effects, turbulence in the reconnection zone, acceleration 
of particles, and comparison to experiments. It has become clear that, espe- 
cially for the acceleration of particles, fully kinetic models and simulations are 
inevitable. Up to now most kinetic simulations make use of the Particle in 
Cell method (PIC), and very few deal with a direct integration of the Vlasov 
equation. In order to validate PIC simulations of collisionless magnetic recon- 
nection, it is necessary to simulate the same problems using different schemes, 
and Vlasov schemes are an important alternative method. 

The second important system currently investigated is the problem of colli- 
sionless shocks. In the framework of MHD-equations these shocks appear as 
singular points, where the macroscopic quantities show a discontinuity. To 
understand the inner structure of shocks, kinetic models are necessary since 
only they can describe the underlying dissipation processes. For some spe- 
cial geometries together with low Mach numbers there are time independent 
analytical solutions to the kinetic equations, corresponding to stable shock 
conditions. For higher Mach numbers these solutions break down, and a time 
dependant behaviour appears. This leads to excitation of wave modes with 
a possibility of particle acceleration. This time dependant behaviour is not 
accessible with analytical methods, and numerical approaches have to be at- 
tempted. 

As mentioned above, in practise there are two different numerical approaches 
to solving the Vlasov equation. In the Particle in Cell (PIC) approach trajecto- 
ries of individual representative particles are followed through the electromag- 
netic fields. The fields are given on a numerical grid, while the positions and 
velocities of the simulation particles can be any value. Due to the finite num- 
ber of particles, the PIC method suffers from considerable numerical noise. A 
related problem is the fact that only those parts of the distribution function 
can be calculated with high precision which contain many particles in a phase 
space volume. In particular the high energy tails of the distribution function 
cannot be resolved. 
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The other numerical approach integrates the Vlasov equation directly on a 
high dimensional numerical grid in phase space. These Vlasov schemes do not 
suffer from any numerical noise. The tails of the distribution function can be 
modelled with high accuracy and deviations from a Maxwellian distribution 
can be pinpointed. These advantages are traded against higher computational 
effort of the Vlasov-codes as compared to PIC-codes. Additionally, care has 
to be taken for choosing the integration scheme of Vlasov's equation, taking 
into account the hyperbolic nature and corresponding integral quantities. 

The distribution function has to fulfil a number of restrictions which are de- 
rived from its physical interpretation as probability density in phase space 
and from the properties of the Vlasov equation. The probability density inter- 
pretation implies that the distribution function has to remain positive at all 
times. One property of Vlasov's equation is that the values of the distribution 
function are transported along the characteristics through phase space with- 
out modification. When starting from a positive distribution function this not 
only implies that the positivity is preserved but also that both upper and lower 
bounds of the distribution function remain unchanged, and that no new max- 
ima or minima are generated. A second property of Vlasov's equation is the 
conservation of phase space density as a consequence of Liouville's equation. 

A numerical scheme cannot satisfy all the above criteria exactly. Therefore a 
number of numerical schemes have been proposed, each of which is a compro- 
mise between different requirements. Spectral codes which solve the Vlasov 
equation in the Fourier domain suffer only from little numerical diffusion, but 
they are mostly limited to periodic boundary conditions, see e.g. [21]. Another 
important drawback of these methods is that they do not preserve positivity, 
let alone the number of extrema in the distribution function. Eulerian solvers, 
on the other hand, allow for non-periodic boundary conditions and can be 
made to preserve positivity and the maximum principle. They are, however, 
slightly more diffusive than spectral codes. A recent comparison of Eulerian 
solvers can be found in [22]. In this work we will use a flux conservative and 
positive scheme [23] which obeys the maximum principle and suffers from 
relatively little numerical diffusion. 

The integration of the Vlasov equation has to be performed simultaneously 
to the evaluation of the Maxwell equations. The dynamics of the full electro- 
dynamic fields imposes an additional criterion for the time step used in the 
simulation. Since electromagnetic waves can travel through the system, the 
time step has to be chosen such that the speed of light c is resolved on the 
numerical grid: At < cAx, where At is the time step and Ax is the grid reso- 
lution. One way to avoid this is to use an electrostatic model. This is, however, 
only applicable in special situations where the self generated magnetic field 
can be neglected. Here we will use the Darwin approximation of Maxwell's 
equations, which follows from a rigorous expansion of the full Maxwell equa- 
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tions in orders of v 2 /c 2 , where v is a characteristic velocity of the system. 
In the framework of the Darwin approximation, the purely electromagnetic 
modes are neglected but electrostatic, magnetostatic, and inductive fields are 
still considered. Darwin's approximation has been widely used with Particle 
in Cell simulations [24,27], however, it has not found its way into Vlasov sim- 
ulations yet. An alternative method is the use of an implicit time stepping as 
it has been implemented by [14] in the Celested3D code. 

The next section will present the basic equations together with their nor- 
malisation. In section 3 we will introduce the Darwin approximation of the 
Maxwell equations. Section 4 will give an overview of the one dimensional 
flux conservative scheme used for integrating Vlasov's equation, while section 
5 will describe the time splitting schemes that provides the generalisation to 
the 5-dimensional phase space. In section 6 results are presented, and section 
7 contains some concluding remarks. 



2 Basic equations 



In this and the next section we want to present the basic equations and ap- 
proximations of our model. The aim is to simulate Vlasov's equation 

^+v.V/ fc +*(E + vxB).V v / t = 0. 

at m k 



Here /^(x, v,t) is the distribution function of species k. In this work only 
singly charged ions and electrons k = i,e are considered, although the code 
allows arbitrary species. The quantities q k and are the charge and the mass 
of the particles of species k. The Lorentz force depends on the electric and 
magnetic fields. These are in general given by Maxwell's equations, 

VxE^-f, (!) 

-VxB = e ^+j, (2) 
fi dt 

V • E = -p , (3) 
V-B = , (4) 

where B is the magnetic and E the electric field. In section 3 we will present 
the Darwin approximation of Maxwell's equations, which is used to solve the 
electromagnetic fields in our simulation code. The charge density p and the 
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current density j are given by the moments of the distribution function, 

P = E* / /fc(x, v)tfV , (5) 
k J 

j = S?*/ v/fe(x,v)d 3 v . (6) 



2. 1 Normalisation 



Here we want to present the normalisation of the Vlasov-Maxwell system 
equations. For this we introduce normalising parameters Aq, where A stands 
for any of the physical quantities. The normalised quantities A are then given 
by A = A/A . 

The un-normalised characteristics of Vlasov's equation for species k are given 
by 

£ = it(E + vxB) . (8) 

at ink 

Here we have used the symbols x and v to denote the characteristics x = 
x(t,xo,vo) and v = v(t, xn,vn). These have to be distinguished from the 
independent variables x and v in Vlasov's equation. In the following the ap- 
propriate meaning should be clear from the context. The form of the above 
equations is not modified by the normalisation, 

dk 

Z = V ' (9) 



S = ifE + vxB) . (10) 



Only the individual charge to mass ratios qk/irik are relevant parameters of 
the system. On the other hand, Maxwell's equations simplify to 

VxE^, (11) 

v*6 = o«(f +3) , <u) 

V-E = p, (13) 
V-B = . (14) 

Here a = vq/c is the ratio of the normalisation velocity over the speed of light. 
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We choose m = rrii and go = e , where is the ion mass and e is the unit 
electron charge. B and n remain free to choose. Then xq = A, is the ion 
inertial length, t = is the inverse ion gyro frequency = ^), and 

Vq = va = B / y/fiQiniUQ is the Alfven velocity. 

We still have the freedom of choosing the magnetic field B . If we choose B 
to be a characteristic magnetic field magnitude in the system, then a gives the 
ratio of Alfven velocity to the speed of light in the system. On the other hand 
we can choose B such that a = 1. Then the local magnetic field magnitude in 
the system determines va/c. In this work we will choose the first, so that the 
field magnitude in the system remains around unity and the free parameter a 
can be used to set the magnetic field strength. In the following we will drop 
the hat-notation (") for the normalised quantities. 



3 Darwin approximation 

The flux conservative integration scheme of the Vlasov equation, which will be 
presented in section 4, is not restricted by a CFL condition on the time step. 
However, when combined with the Maxwell equations for the electromagnetic 
fields a CFL condition is introduced by the time integration of the fields on 
the numerical grid. This implies that the fastest electromagnetic wave mode, 
i.e. the vacuum mode, has to be resolved on the grid, 

At < ^ , (15) 

where c is the speed of light. This condition imposes severe restrictions on the 
time step. In many applications the electromagnetic vacuum modes are not 
important. The standard solution to this problem is the electrostatic approx- 
imation. In this approximation only Poisson's equation needs to be solved for 
the electric field. A magnetic vacuum field can also be superimposed. The in- 
fluence of the plasma on the magnetic field is, however, completely neglected. 
To account for this influence, and thus the possibility of magnetosonic wave 
modes, the Darwin approximation is commonly used in particle simulations. 
This approximation can be derived from the Maxwell equations in an ex- 
pansion in orders of v 2 /c 2 , where v is some characteristic velocity. Assuming 
v 2 <C c 2 this leads to a set of equations where the vacuum modes are elimi- 
nated but all other wave modes are retained. Darwin's approximation can be 
used when the velocities are small compared to the speed of light and there 
is no energy transported by the electromagnetic radiation. Some applications 
are the study of magnetic reconnection, e.g. in the earths magnetotail [25], 
or the investigation of high intensity charged particle beams [26]. Naturally 
it is not capable of describing phenomena where the electromagnetic vacuum 
modes play a major role, e.g. in laser-plasma interaction. 
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Darwin's approximation starts from a separation of the electric field into a 
longitudinal and a transverse part [27] 



E — El + , 



(16) 



with 



V x E L = and V • E T = . 



(17) 



The normalised Maxwell's equations are then approximated by 



V-E L 
VB 

VxE T 
V x B 



P , 


-d t B , 

a 2 (d t E L +j) 



(18) 
(19) 
(20) 
(21) 



The approximation only appears in eq (21), where only the longitudinal part 
of the displacement current is taken. All other equations remain unchanged. 
In vacuum, equations (20) and (21) are now decoupled and no purely electro- 
magnetic modes can appear. 

An advantage of using Darwin's approximation instead of the full Maxwell 
equations is the fact, that the equations can be solved without performing a 
time integration step. All the electromagnetic fields can be calculated from 
the moments of the distribution at a given time together with the boundary 
conditions for the fields. This will be presented in the following. 

Poisson's equation for the electrostatic potential 



immediately gives the longitudinal electric field. Taking the curl of (21), to- 
gether with (17) and (19, results in 



giving three Poisson equations for the components of the magnetic field. To 
calculate E T the curl of eq (20) is taken and eqs (21) and (17) substituted, 
giving 



To eliminate the time derivative of the current density, j is expressed as mo- 
ment of the distribution function, 



A$ = -p with E L = -V$ 



AB = -a 2 V x j , 



AE T = dfV x B 




(22) 
(23) 




k 
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Here the index k sums over all particle species. Now, Vlasov's equations for 
the different species are substituted, 



^j = -£Vp fc <vv) fc + £^E + £^<v B 

k k m k k m k 

Here the pointed brackets (.)& denote averaging with the distribution function 
fk- The electric field appearing on the right hand side is, of course, comprised 
of the longitudinal and the transverse part. While the longitudinal component 
is already known from the charge density, it is the transverse component that 
is being calculated here. Inserting <9J into eq (23) and introducing the local 
plasma frequency 



2 _^q_kPk ? (24) 



this leads to 



AE T - « 2 cj 2 E t = 

-EVp,(vv) fc + 5:^E L + 5:^(v) fc xB-V(« 2 ^0) . (25) 

k k m k k m k 

This equation is a Helmholtz-equation for each component of Ey. The last 
term on the right hand side can be deduced from the condition V • E r = 0, 
since it only adds a curl free component to Ey. It can thus be omitted to 
calculate E T , with 

AE r - aWV T = -J2 V Pfc (vv) fe + £ + £ ^(v) fe x B . 

k k m k k m k 

Et is then projected onto its divergence free part by calculating with AG = 
V • E T to give E T = E T — V6. 

Altogether this sums up to 8 elliptic equations. However, the time spent in 
solving these equations, compared to the integration of the 5 or 6 dimensional 
problem, is negligible. This situation is different in PIC simulations. One has 
to compare about 30 particles per cell used in standard PIC simulations to 10 3 
- 30 3 mesh points for the resolution of the velocity space used in Vlasov codes. 
For this reason the computational effort of solving the above elliptic equations 
has an effect on the total computational time in standard PIC simulations, 
but not in Vlasov simulations. 



4 Flux Conservative Scheme 



In this section we briefly present the numerical scheme used for integrating 
the Vlasov equation. This scheme has originally been presented in Filbet et 
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al. [23]. The scheme uses a flux conservative formulation and is based on a 
third order reconstruction of the primitive of the distribution function using 
a fixed stencil. 

To formulate the scheme we start from the observation that in Hamiltonian 
systems the values of the distribution function are transported along the char- 
acteristics 

f(U)^f(X(s,t,0,s). (26) 

Here X(s,t,£) denotes the characteristic with parameter s that satisfies 
X(t,t,£) = £, where in general £ = (x, v). For the rest of this section, we 
restrict the calculations to one dimension. The generalisation to the higher 
dimensional system will be given later. We thus assume ( 6i Then we can 
integrate (26) and obtain a propagator from time t n to time t n+1 

X(t n ,t n +\x l+1/2 ) 

J f(x,t n+l )dx = I f(x,t n )dx. (27) 



Here Xj_i/ 2 and x i+ i/ 2 are the boundaries of the numerical grid cell i. The dis- 
cretisation of the distribution function is now suggested by the above equation. 
The values /" on the numerical grid represent the cell integrals 

K = -t f f(xX)dx. 

We also define the flux through a cell boundary at Xj+i/2 during the time 
interval [t n ;t n+1 ] 

= / f(x,ndx 

X(t»,t»+i,x i+1/2 ) 

Then equation (27) can be expressed in the following form 

= *?-i/2 + f?- ^i/2- (28) 

This equation still follows exactly from Vlasov's equation. However, the flux 
$ has to be calculated, and this requires an approximation of the distribution 
function. 

In the framework of the third order positive and flux conservative scheme, as 
presented by [23], the primitive of the distribution function is approximated 
using a four point stencil. Let 

F(x,t n )= f f(x',t n )dx', (29) 

Jxr, 
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then it follows exactly that 



F(x t+l/2 ,t n ) = AxJ2fJ! = F l n . 

k=0 



(30) 



The primitive is given exactly on the cell boundaries. To approximate F(x, t n ) 
in the cell interval [xj_i/ 2 ; x i+1/2 ], the four points {2^-3/2, Xi-1/2, ^+1/2, ^+3/2} 
are used. Taking the derivative of the primitive, we recover the approximation 
of the distribution function 



Ux) = fi + 



2(x- Xi) (x - Xi_3/ 2 ) 



6Ax 2 

+ (x- Xi-!/ 2 ) (x - Xi+1/2)] (fi+1 - fi) 

~ 6A^I 2 ( x ~ Xi) (^-^+3/2) 

+ (a: - Xi-1/2) (a: - Si+1/2)] (/i - fi-i) ■ 



(31) 



Here limiters ef are introduced. The limiters are chosen in such a way as to 
limit the values of the distribution function to a fixed interval < f&(x) < /oo, 
where /oo is the maximum of all fi. The ef can be written as 



,+ 



min(l;2/,/(/ m -/ t )) 



if fi+i -fi>0 



and 



min(l; 2(/ 0O - /,)/(/, - f i+1 )) if - /< < 

min(l; 2(/ 00 - /,)/(/, - /*-i)) if /, - £-1 > 
min(l;2/ i /(/ i _i-/ i )) if fi - fi-i < ' 



(32) 



(33) 



If we determine j such that £7-1/2 < -^(^ n , t n+1 , 2^+1/2) < 2^+1/2, we get 
$ ,+i/ 2 = (1 - 5) [fj + (l/QW + l)e + (/,+i - fi) 



max(i,j)-l (34) 

(1/6)5(5 -2)e-(fj -f j - 1 ))]+ Ax A, 

fc=min(i,j')+l 



where 5 = X(t n , t n+1 , x i+ i/ 2 ) - £7-1/2- 



5 Time splitting and integration of characteristics 



The above scheme was presented in one dimension. To perform two dimen- 
sional Vlasov-simulations including magnetic fields, the scheme has to be gen- 
eralised to a total of 5 dimensions. Califano et al. [28] proposed a scheme in 
which the individual steps for the different dimension are carried out using a 
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second order time splitting scheme. The scheme is an extension of the time 
splitting scheme presented by Cheng and Knorr [29]. For the two spatial di- 
rections x and y, the projection of the characteristics onto the axes is simply 
evaluated from by the corresponding velocity components v x and v y . In this 
way the two space dimensions are independent of each other and can be per- 
formed sequentially, 

T* = T x T y (35) 
where T\. denotes a time-integration step in one dimension with the direction 
k. In contrast to this, the velocity components are not independent of each 
other due to the magnetic field. To create a second order scheme Califano et 
al. [28] formulated a straightforward second order time splitting scheme, 

At At At At At At 

T v (At) = T vx {—)T vy {—)T vx {—)T vz {At)T vx (—)T vy {—)T vx (—). (36) 

This propagator in the velocity directions was then combined with the prop- 
agator in the space directions to create the full scheme, 

T full (At) = T x (At/2)T v (At)T x (At/2). (37) 

In the following we will refer to this method as the time splitting method. 

We implemented and tested this scheme, but found it to suffer from substantial 
inaccuracies, as shown in section 6.1. Therefore, a new scheme for the inte- 
gration is proposed in this work which we call the back substitution method. 
Here, the integration of Vlasov's equation in the three velocity dimensions is 
still split into three separate steps, one for each direction. The difference to 
the above scheme is the way the integration of the characteristics is separated 
out into the three substeps. This implies that the concrete implementation 
of the scheme proposed here depends on the method of integration of the 
characteristics. We will therefore quickly present Boris' scheme for integration 
of characteristics in a magnetic field [27,30,31]. Boris' scheme is widely used 
in Particle codes and has the advantage that it is a second order integration 
scheme in which the magnetic field does not cause a change in the kinetic en- 
ergy. The integration step is formulated as an implicit finite difference scheme 

n+l _ n / n+1 , n \ 

- x B. (38) 



At m\ 2 

The electric and magnetic forces are separated, 



v" + ^E, (39) 
2 m 



v v" 1 — — E, (10) 

2m' y 1 

leading to 



v+ — 



At 2m 
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The transformation from v to v + is a pure rotation with an angle 6, where 



9 AtqB 



(42) 




To implement this rotation the vectors t and s are defined 



t 



At qB 

2 m ' 



s = 



1 + t 2 ' 



2t 



(43) 



Then the rotation is performed in two steps 



v' = v + v x t 



(44) 



and 



v + = v + v' x s. 



(45) 



This scheme supplies v™ +1 = (v v v ™ +1 ) in terms of v n = (t>™, t>™, t>"). 

For the integration of Vlasov's equation in three dimensional velocity space the 
three individual integration steps are performed in turn, thus, the integration 
of the characteristics also has to be split into three separate steps. For this, 
two things have to be considered. Firstly, when integrating one component 
of the velocity, we take into account the shifts of the distribution function 
that have already been performed. Secondly, we then have to observe the 
order of integration, which is reversed with respect to the integration of a 
particle trajectory. This follows from the fact that, for the Vlasov scheme, the 
characteristics have to be traced backwards in time. 

For clarity, we will present the back substation method for forward integra- 
tion of the characteristics. Later, we will reverse the order of the substeps for 
the backward integration needed in the Vlasov scheme. For forward integra- 
tion, first, the integration in v x is performed, as described in the discussion 
of Boris' scheme above. For the integration of the characteristics this means 
v x +1 = v x +1 ( v xi v y-> v z )• When the integration along the i^-direction is per- 
formed, in the next step, the shift along the i^-direction has already been 
performed. For the integration of the characteristics this means that a scheme 
calculating v r y l+1 = i>" +1 (i>£ +1 , i>™, i>") is needed. In the above notation, it is 
sufficient to reformulate the integration to give Vy = Vy(v^,v~,v~). This can 
be achieved using simple algebraic manipulations and eliminating v~. For the 
last integration step along the v 2 -direction, one can find in the same manner 
a scheme giving t>+ = Vy,v~). The details of this scheme are presented 

in appendix A. The scheme has also been implemented and it was found to 
be less diffusive than the second order time splitting scheme (see section 6.1). 
As stated above, the characteristics have to be traced backwards in time. To 
this end, we reverse the order of integration, without changing the above for- 
mulas; i.e. we first shift the distribution function in the t> 2 -direction using 
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v t = v t(, v xi v yi v z )? then in the f y -direction using v+ = Vy(v. 
finally in the i^-direction using t>™ +1 = t>" +1 (t> ™, v ™, v "). 



,v y ,v z ), and 



The splitting of the integration into individual steps introduces another prob- 
lem, independent of the splitting scheme used. The limiters ef and guar- 
antee only that the distribution function is positive and limited by from 
above for incompressible transport equations. The Vlasov equation together 
with any set of equations for the electromagnetic fields presents such a sys- 
tem. Due to the separation of the integration along the different directions 
in velocity space together with the forces originating from the magnetic field, 
the single integration steps are, however, not incompressible. This holds for 
both the splitting and the back substitution scheme. The first direction of 
integration might compress the distribution function in the v x direction, while 
the successive integrations decompress the distribution function in the v y and 
v z directions. In the end the total compression will always vanish. The max- 
imum value of the distribution function might however increase during an 
intermediate step. 

We present two types of calculations. In the first we omit the limiter from 
above completely, allowing the distribution function to rise uncontrollably. 
The limiters ef and e~ are then given by 



This means, the maximum principle, that was fulfilled by the original scheme, 
is now no longer satisfied. In the second type of simulation we keep the original 
limiters (32) and (33) but determine foo to be the global maximum of the 
distribution function after each intermediate integration step for each direction 
in velocity space. 



6 Numerical tests 

The third order positive flux conservative method has been extensively tested 
in one dimension and in the electrostatic limit [23]. Here, we first want to 
present results of the integration scheme in three dimensional velocity space 
with a given magnetic field. We compare the results of the two integration 
schemes and a standard second order finite difference scheme. The latter is 
included in the comparison since it has been used for Vlasov simulations of 
reconnection by various authors. 




(46) 



and 




if fi ~ fi-i > 
fi)) if /i - < 



(47) 
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In the second part of this section we present results of reconnection simula- 
tions. These have been carried out using both the flux conservative scheme 
with the back substitution method and the finite difference scheme. Results 
of the two are compared. 



6.1 Gyro motion 



To test the quality of the Vlasov integration scheme together with the inte- 
gration of the characteristics a simple test system was simulated. In this test 
only one positively charged species was simulated in a constant magnetic field 
B = B e z , with B — 1. The initial distribution function was taken to be a 
shifted Maxwellian 

/(x,v) = exp[-(v-v ) 2 ], (48) 

where v = v oe x is a constant velocity in the x-direction. There was no spatial 
variation, and so a simple gyro motion of the thermal peak in velocity space 
is expected. However, due to numerical errors this peak will change its shape. 
These numerical errors only appear in the presence of a magnetic field. In the 
pure one dimensional advection problem only small errors compared with the 
exact solution are found. 

We have simulated the gyro motion of a thermal peak using the two different 
integration methods for the characteristics, different time steps, and vary- 
ing resolutions of the grid in the velocity dimensions. Additional simulations 
have been carried out using a simple finite difference scheme. The scheme is 
obtained from eq (2) by substituting all partial derivatives by their centred 
finite difference approximation over one grid cell. A third order Runge-Kutta 
method is used for the time step. In Figure 1 the v x -v y distribution function is 
shown for the case of N v = 30 and At = 27rf2 c /200 using the flux conservative 
scheme with the back substitution algorithm for integration of the character- 
istics. The z-scaling is arbitrary but identical in the three sub-plots of Figure 
1. Initially (a), the distribution function is set to the shifted Maxwellian. At 
t = 2nfl c (b), the peak has performed one full gyration. A slight dissipa- 
tion can be observed. The dissipation has increased after five gyration periods 
t = 107rf2 c (c). With the above time step, this correspond to a total of 1000 
time steps. 

The dissipation in phase space can be quantified by taking the second moment 
of the distribution function ((v — (v)) 2 } = (v 2 ) — (v) 2 which represents the 
thermal energy in the system. These have been calculated over a time of 
five gyro-periods for different simulation time steps, different grid resolutions 
of the phase space, and different time splitting algorithms. In Figure 2 the 
thermal energies are shown for both the back-substitution algorithm (upper 
panel) and for the time splitting algorithm (lower panel), both taken without 
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Fig. 1. Gyration of a Maxwellian peak using the back-substitution algorithm with 
N v = 30 and At = 27rfi c /200 (a) initially, (b) after one gyration period and (c) after 
five gyration periods. The scale on the z-axis is in arbitrary units but identical in 
the three diagrams. 

limiting the distribution function from above. For a small time step At = 
0.01f2 c /7r, i.e. a time-resolution of 200 time steps per gyration, calculations 
using N v — 10 and N v = 30 grid points per velocity dimension have been 
carried out. In the case of N v = 10 and using the back-substitution algorithm, 
the thermal energy rises more than a factor of 4.5 during the simulated time 
period. The amount of diffusion can be reduced considerably by increasing the 
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Fig. 2. The thermal energy of the distribution function (v 2 ) — (v) 2 during five 
gyro-periods for different time-steps and different phase space resolutions. The 
upper panel shows the curves resulting from the back-substitution algorithm, the 
lower panel those resulting from the time splitting algorithm. No limiting from above 
was applied 

grid resolution. At N v = 30, the increase in thermal energy reduces to around 
20% over the simulated five gyro-periods. 

For N v = 10, an increase in the time step by a factor of 10 does not modify the 
results by a great amount. The already large diffusion values remain similar 
although somewhat lower. In contrast, for the N v = 30 case, the thermal 
energy at the end of the simulation is risen only by less than 5% from the 
starting value, for a time step of At = 0.1Q c /tt, i.e. 20 time steps per gyration. 



16 



This means that an increase in the time steps leads to better results for the 
thermal energy. An increase in grid resolution further improves the result and 
reduces the error to less than 2%. 

The results using the time splitting scheme (lower panel of Figure 2) differ 
considerably from the previous results. Again, the thermal energy in the small 
time step case with N v = 10 rises strongly. This now gets substantially less 
accurate by increasing the time step to 20 steps per gyration. Here, the thermal 
energy increases by more than a factor of 12! The best result is achieved by 
taking a small time step and a grid resolution of N v = 30. Instead of a slight 
increase in thermal energy, a slight decrease can be observed here. Increasing 
the time step to 20 steps per gyration degrades the results again. A strong 
increase in the thermal energy is observed, which is even worse for higher grid 
resolutions. 

In Figure 3 the same results are shown again, but this time with the lim- 
iting from above switched on. These results differ only marginally from 
those obtained without the limiter in place. Note, that the limiter ensur- 
ing positivity of the distribution function is always in place. Summarising, 
one finds that the back-substitution method gives superior results over the 
time-splitting method. Furthermore in terms of thermal energy, the back- 
substitution method rewards larger time steps (i.e. less computational time) 
with higher accuracy. 

Instead of looking at the thermal energy, which is related to dissipation in 
phase space, one can also investigate the gyro-motion of the Maxwellian peak. 
From the distribution function the averages (v x ) and (v y ) are taken and plotted 
against each other through time. In the exact case, this should result in per- 
fectly circular motion with exactly one turn during a gyration period 2n/Q C - 
The numerical Vlasov scheme does not, however, produce exact results and 
thus one will observe some form of spiral. An example of this spiral is plotted 
in Figure 4. 



After some time one can observe errors in both the absolute magnitude of 
the velocity and in the phase of the gyro-motion. For the same parameters 
as above, the magnitude of the velocity after five gyro-periods is given in 
table 1, for the various schemes and numerical parameters. Table 2 shows the 
corresponding phase errors in radians. 

The back-substitution shows better results both for the large and the small 
time steps than the time splitting method. The best overall result is achieved 
by the back-substitution method in conjunction with a small time step and a 
grid resolution of N v = 30. Only a marginal phase error can be observed. There 
is almost no error in the average velocity. A decrease of the grid resolution 
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Fig. 3. The thermal energy of the distribution function (v 2 ) — (v) 2 during five 
gyro-periods for different time-steps and different phase space resolutions. The 
upper panel shows the curves resulting from the back-substitution algorithm, the 
lower panel those resulting from the time splitting algorithm. Limiting from above 
was applied 

results in a positive phase error and a decrease of the velocity magnitude error. 
In contrast, an increase of the time step causes a rise in the velocity magnitude 
together with an increase of the phase error. Here now, an increase in the grid 
resolution from 30 to 60 does not improve the results when using a large time 
step. For a resolution of 10 and a large time step only the phase error increases. 
One can note, that in all these cases the phase error is negative. This means 
that the simulated gyro motion always lags behind the physical gyro motion. 
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Fig. 4. The gyro-motion of the Maxwellian peak expressed in (v y ) against (v x ) 
over five gyro-periods. After five gyro-periods the mean velocity has an error in 
magnitude and phase from the starting value. 
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Table 1 

Magnitude of the Maxwellian peak after five gyro-periods for different time-steps, 
different phase space resolutions and different integration methods. The magnitudes 
are normalised to the initial height of the peak. 



The errors in phase and velocity magnitude are generally larger when using the 
time splitting method. With the large time step the magnitude of the velocity 
goes up by approximately 2 to 4 times the initial velocity. The velocity error 
is lower when using the small time step; however, for the high resolution the 
error is still 20%. Only the low resolution, large time step case has a relatively 
small error in the magnitude. This is accompanied by a large error in phase. 
For resolutions 30 and 60 the phase error is better than the one found with 
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Table 2 



Phase error of the Maxwellian peak after five gyro-periods for different time-steps, 
different phase space resolutions and different integration methods. 

the back-substitution method. The error in the velocity magnitude here is on 
the other hand, unacceptable. In both cases one finds that switching on the 
limiter does not modify any of these results considerably. 

Finally, comparing this with the results from the finite difference scheme one 
finds that the flux conservative approach is far superior in all cases. Although 
the finite difference schemes show small amplitude errors, phase errors are se- 
rious. In addition, for the large time steps the finite difference scheme becomes 
completely unstable (in agreement with the CFL condition). 

6.2 Reconnection 

A test of the complete Vlasov-Darwin system has been performed on a mag- 
netic reconnection setup, including full ion and electron dynamics. As pointed 
out in the introduction, research on magnetic reconnection is still one the most 
challenging topics in collisionless plasmas. Kinetic simulations using a Vlasov 
code have been carried out, for example by Silin and Biichner [18]. Since these 
simulations were carried out using a finite difference scheme, we will here again 
present the results of both finite difference and the flux conservative scheme. 
We used a simulation box consisting of 50x100 grid cells in space with a grid 
spacing of Ax = 0.1. The total size of the simulation box is L x = 10 by 
L y = 5. There are 10x10x10 grid cells in velocity space. Although the resolu- 
tion in velocity space is — as shown above — insufficient for the gyro motion, 
especially for the finite difference scheme, this low resolution is chosen for a 
number of reasons. Firstly, it keeps the computational effort to a reasonable 
magnitude. Secondly, this corresponds to the phase space resolution chosen 
by Wiegelmann and Biichner [19] and thus allows some comparison. Thirdly, 
we think that the pure gyro motion test is a very strict test for the numerical 
scheme. In the most regions of the reconnection simulation, the magnetic field 
is balanced by some electric field or by the density gradient in such a way that 
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Fig. 5. The reconnected magnetic flux $ over time. The solid line is the result of 
the second order finite difference scheme. The dashed line is the result of a fourth 
order finite difference scheme. The dashed-dotted line is the result from the flux 
conservative scheme. 

the gyro motion does not occur the way as presented in the previous section. 
The mass ratio between electrons and ions has been chosen to mj/m e = 16 
and the time step is 

A( = 2 ' 5 ■ 10-3 = 25b (49) 

where Q ce is the electron Larmor frequency in the unit magnetic field. Simu- 
lations with half this time step have been performed as convergence test (not 
shown here) but no substantial deviations have been found. 

The initial conditions are chosen as two opposite current sheets. Each cur- 
rent layer has the well known Harris sheet profile [1] together with a small 
perturbation. The distribution function is given by a shifted Maxwellian 

fi,e{x,y) = p(x,y)exp (yl + vl + (v z -v^) 2 )^j , (50) 

with Vq = ±1. The particle density p(x,y) is given by 

V) = —T7 ZVTTY 1 + £ cos(^- . (51) 



cosh((y-y^) /A) 



The two Harris sheets, positioned at y + = L y /A and y = 3L y /A, carry op- 
posite current, and the perturbations have a relative shift along the x-axis 
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Fig. 6. The out of plane magnetic field component B z in the simulation box at 
t = 4.1 using the back-substitution method (upper panel) and at t = 3.1 using the 
finite difference scheme in 2nd order (middle panel) and 4th order (lower panel). 

of tt/L x , given by x + = and x~ = 1/2. The perturbation has a value of 
e = 0.05. Periodic boundary conditions are used both in the x and y direc- 
tion. The X-points will then develop at X + = (5, 1.25) and X~ = (0, 3.75). 

Figure 5 shows the temporal evolution of the reconnected flux. The curves 
are quite similar keeping in mind that the resolution is marginal. The major 
difference is the slightly earlier onset of reconnection for the finite difference 
runs. 



If one wants to compare the schemes, one has to choose the time for the dif- 
ferent simulations separately such that the reconnected flux is the same for all 
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Fig. 7. The magnitude of the Hall term |(j x B) z | in the simulation box at t = 4.1 
using the back-substitution method (upper panel) and at t = 3.1 using the finite 
difference scheme in 2nd order (middle panel) and 4th order (lower panel). 

simulations. Therefore, we choose a time at the beginning of the reconnection 
phase which corresponds to t — 4.1 for the flux conservative method and to 
t = 3.1 for the finite difference schemes. 

Figure 6 shows the perpendicular magnetic field B z . The upper panel shows the 
result of the flux conservative simulation while the mid and lower panel display 
the results of the finite difference schemes (mid: second order, lower: fourth 
order). One can clearly see the quadrupolar structure of the magnetic field 
component at the X-point which has been observed in previous simulations of 
the reconnection process. There is a slight difference in the distribution of the 
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magnetic field component between the simulations. One significant difference 
is the magnitude of B z . The values are considerably larger in the second order 
finite difference scheme than in the flux conservative scheme and fourth order 
finite difference scheme. 

In Figure 7 the magnitude of the ^-component of the Hall-term |(j x B) 2 | 
is shown for the different schemes. This quantity can be taken as a measure 
for the importance of the Hall term when comparing to MHD models. All 
simulations show qualitatively a similar behaviour. From the results of the 
flux conservative simulation one can see that the z-component of the Hall 
term is largest outside of the X-points.l 



7 Conclusions 

A 5 dimensional Vlasov code using the Darwin approximation of Maxwell's 
equations has been presented. While in the past a lot of development has gone 
into developing numerical schemes for one dimensional electrostatic Vlasov 
codes, there is a growing need for kinetic codes that include the effects of the 
magnetic field. This is needed for both the simulation of astrophysical as well 
as laboratory plasmas. For a large number of problems the electromagnetic 
vacuum modes contained in the full set of Maxwell equations is not of major 
importance for the physical processes. These, however, pose severe restric- 
tions on the simulation parameters, due to the CFL criterion. The Darwin 
approximation resolves this problem by eliminating the vacuum modes with- 
out reverting to the electrostatic limit. In contrast to other approximation of 
Maxwell's equations, the Darwin approximation can be shown to consistently 
emerge from an expansion of the full equations in orders of Vq/c 2 . The Darwin 
approximation has been used extensively in particle in cell simulations but has 
not found it's way into Vlasov codes until now. Within the framework of the 
Darwin approximation the purely electromagnetic modes are cancelled, but 
the modes that rely on the plasma reaction to the magnetic field are retained, 
such as the magnetosonic or Alfven modes. 

For the integration of the Vlasov equation in time a recently developed flux 
conservative scheme [23] has been used. The scheme, which was originally pro- 
posed for a one dimensional electrostatic system had, to be generalised for the 
higher dimensional phase space required for the Vlasov-Darwin system. The 
generalisation of the one dimensional system is ambiguous and two different 
schemes have been compared: the straightforward time splitting scheme [28] 
using a generalised time splitting iteration was shown to suffer from substan- 
tial inaccuracies; the back substitution scheme, proposed in this work, gives 
far better results not only in terms of an accurate reproduction of the gyro 
motion but also concerning the controllability of the errors. Results were also 
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compared to finite difference simulations which were shown to be greatly in- 
accurate in comparison to the flux conservative scheme. 

Finally, a simulation of magnetic reconnection has been performed to test the 
full Vlasov-Darwin system. For comparison, the same system has also been 
simulated using the finite difference scheme. The finite difference scheme has 
been chosen as a reference scheme since this is the most common approach for 
simulation magnetic reconnection using a Vlasov code. Substantial improve- 
ments were achieved using the code presented here. The results especially of 
derived quantities, did not suffer from numerical fine scale distortions. It was 
shown that the results of the scheme presented here can be treated as much 
more reliable that those of a finite difference scheme. 

The code presented here may be used for the investigation of magnetic recon- 
nection processes and for the simulation of nonrelativistic collisionless shocks. 
The main advantages, as compared to PIC simulations, is the absence of nu- 
merical noise and the access to the full distribution function, especially the 
high energy tails. In PIC simulations these high energy tails generally suffer 
from bad statistics. For this reason it might be of interest to use Vlasov codes 
to simulate particle acceleration. 

One should note however that, because of the large computational time, the 
applicability of Vlasov codes at this time stays restricted to two dimensional 
models. Again, in comparison to PIC codes, Vlasov codes are slower by a fac- 
tor of a hundred. Despite this comparatively large numerical effort, Vlasov 
simulations are still valuable since they complement other simulation tech- 
niques. 
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A Back-Substitution Method for Integration of Characteristics 

Here we present the formulae for the integration of the characteristics in the 
velocity space using the back-substitution method. The three integrations in 
v x , v y and t> 2 -direction are carried out individually. The underlying equations 
have been presented in section 5 but are in an implicit form for our purpose. 
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Here we want to give explicit formulae for the three steps 



< +1 = < +1 «,<,<), (A.l) 

< +1 = < +1 (< +1 ,<,<), (A.2) 

< +1 = < +1 « + \< + \<). (A.3) 

Since the bijections between v n and v~ on one hand and v n+1 and v + on the 
other hand are trivial (see eqs (39) and (40)) it is sufficient to formulate the 
three steps 

vt = v+(v-,v-,v~), (A A) 

Vy =v+(vZ,Vy,v~), (A.5) 

vt = vt(v2,v+,v~). (A.6) 



By virtue of equations (44) and the x-component of (45), v x is already given 
in the above form. The integration in v x is identical to the integration used in 
the time splitting algorithm. 

In the second integration step Vy has to be determined from (u+ ,v~ ,v~). 
The y-component of equation (45) gives Vy = Vy(v~, v~, v~). However, the 
re-component of that equation can now be solved for v~ to give 

=v x -(vZ,v y -,v;) = ^- (u+ -v-(s z + s y t x ) -v-{s y - s z t x )) , (A.7) 



x 



with 

A x = 1 Syty s z t z (A. 8) 

With this we get 

vj(v2,v-,v-) =vf(v-(vZ,v-,v-),v-,v-). (A.9) 

In the same manner the z-component of equation (45) gives = 
v+(v~, v~, v~). Then the x and the y-components of that equation are used 
to solve for v~ and v~ 

x y 



v„ = 



(1 - S x t x - S Z t z ) + (Sy - S Z t X )V z )-(S z + Syt X ) - (S x + S Z ty)V z ) 

(1 Syty s z^z)(l s x^x s z^z) ~l~ ( s z Syt x ){s z S x ty) 

(A.10) 

(S Z - S X ty) (U+ + (Sy ~ S Z t X )V~) + (1 - Syty ~ S g t g) ~ (S X + S z ty)V~} 



V = 

(1 Syty S 2; t 2 )(l S x t x S z t z ) + (s 2 + Syt x )(s z S x ty) 



(A.11) 
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Inserting this into v+(v x , v y , v z ) then provides the expression for t>+ = 



References 



[I] E. G. Harris, On a plasma sheath separating regions of oppositely directed 
magnetic field, II Nuovo Cimento 23 (1962) 115. 

[2] D. Biskamp, E. Schwartz, and J. F. Drake, Two fluid theory of collisionless 
magnetic reconnection, Phys. Plasmas 4 (1997) 1002. 

[3] R. F. Lottermoser and M. Scholer, Undriven magnetic reconnection in 
magnetohydrodynamics and Hall Magnetohydro dynamics, J. Geophys. Res. 
102 (1997) 4875. 

[4] M. A. Shay and J. F. Drake, The role of electron dissipation on the rate of 
collisionless magnetic reconnection, Geophys. Res. Lett. 25 (1998) 3759. 

[5] M. A. Shay, J. F. Drake, B. N. Rogers, and R. E. Denton, Alfvenic collisionless 
magnetic reconnection and the Hall term, J. Geophys. Res. 106 (2001) 3759. 

[6] X. Wang, A. Bhattacharjee, and Z. W. Ma, Scaling of Collisionless Forced 
Reconnection, Phys. Rev. Lett. 87 (2001) 265003. 

[7] J. D. Huba and L. I. Rudakov, Hall Magnetic Reconnection Rate, Phys. Rev. 
Lett. 93 (2004) 175003. 

[8] J. Biichner and J. -P. Kuska, Sausage mode instability of thin current sheets as 
a cause of magnetospheric substorms, Ann. Geophys. 17 (1999) 64. 

[9] R. Horiuchi and T. Sato, Three-dimensional particle simulation of plasma 
instability and collisionless reconnection in a current sheet, Phys. Plasmas 6 
(1999) 4565. 

[10] P. L. Pritchett, Geospace Environment Modeling magnetic reconnection 
challenge: Simulation with a full particle electromagnetic code, J. Geophys. 
Res. 106 (2001) 3783. 

[II] P. L. Pritchett, Collisionless magnetic reconnection in a three-dimensional open 
system, J. Geophys. Res. 106 (2001) 25961. 

[12] M. Hesse, M. Kuznetsova, and J. Birn, Particle-in-cell simulations of three- 
dimensional collisionless magnetic reconnection, J. Geophys. Res. 106 (2001) 
29831. 

[13] R. Sydora, Nonlinear dynamics of small-scale magnetic islands in high 
temperature plasmas, Phys. Plasmas 8 (2001) 1929. 

[14] P. Ricci, G. Lapenta, and J. U. Brackbill A Simplified Implicit Maxwell Solver, 
J. Comp. Phys. 183 (2002) 117. 



27 



[15] P. Ricci, G. Lapenta, and J. U. Brackbill, GEM reconnection challenge: Implicit 
kinetic simulations with the physical mass ratio, Geophys. Res. Lett. 29 (2002) 
2088. 

[16] A. Zeiler, D. Biskamp, J. F. Drake, B. N. Rogers, M. A. Shay, and M. Scholer, 
Three-dimensional particle simulations of collisionless magnetic reconnection, 
J. Geophys. Res. 107 (2002) 1230. 

[17] M. Scholer, I. Sidorenko, C. H. Jaroschek, R. A. Treumann, and A. Zeiler Onset 
of collisionless magnetic reconnection in thin current sheets: Three-dimensional 
particle simulations, Phys. Plasmas 10 (2003) 3521. 

[18] I. Silin, J. Biichner, Kinetic instabilities of thin current sheets: Results of two- 
and-one-half-dimensional Vlasov code simulations, Phys. Plasmas 10 (2003) 
1299. 

[19] T. Wiegelmann, J. Biichner, Evolution of magnetic helicity in the course of 
kinetic magnetic reconnection, Nonlinear Processes in Geophysics 8 (2001) 127. 

[20] J. Birn, J. F. Drake, M. A. Shay, B. N. Rogers, R. E. Denton, M. Hesse, M. 
Kuznetsova, Z. W. Ma, A. Bhattacharjee, A. Otto, P. L. Pritchett, Geospace 
Environmental Modeling (GEM) magnetic reconnection challenge, J. Geophys. 
Res. 106 (2001) 3715. 

[21] T. P. Armstrong, R. Harding, G. Knorr, D. Montgomery, Solution of Vlasov's 
equation by transform methods, Advances in Computational Physics 9 (1976) 
29. 

[22] T. Arber, R. G. L. Vann, A critical comparison of Eulerian grid based Vlasov 
solvers, J. Comp. Phys. 180 (2002) 339. 

[23] F. Filbet, E. Sonnendriicker, P. Bertrand, Conservative numerical schemes for 
the Vlasov equation, J. Comp. Phys. 172 (2001) 166. 

[24] C. W. Nielson, H. R. Lewis, Particle-code models in the nonradiative limit, 
Meth. Comput. Phys. 16 (1976) 367. 

[25] P. L. Pritchett, F. V. Coroniti, R. Pellat, H. Karimadbadi Collisionless 
reconnection in two-dimensional magnetotail equilibria J. Geophys. Res. 96 
(1991) 11523. 

[26] W. W. L. Lee, E. Startsev, Qin Hong, R. C. Davidson Electromagnetic (Darwin) 
model for three-dimensional perturbative particle simulation of high intensity 
beams. PACS2001, Proceedings of the 2001 Particle Accelerator Conference, 
IEEE Part vol. 3 (2001) 1906 

[27] C. K. Birdsall, A. B. Langdon, Plasma Physics via Computer Simulation, 
McGraw-Hill, New York, 1985. 

[28] F. Califano, A. Mangeney, C. Cavazzoni, P. Travnicek, A numerical scheme 
for the integration of the Vlasov-Maxwell system of equations, in: Science and 
Supercomputing at CINECA, 2001, p. 456. 



28 



[29] C. Z. Cheng, G. Knorr, The integration of the Vlasov equation in configuration 
space, J. Comp. Phys. 22 (1976) 330. 

[30] O. Buneman, Time reversible difference procedures, J. Comp. Phys. 1 (1967) 
517. 

[31] J. P. Boris, Relativistic plasma simulation-optimization of a hybrid code, in: 
Proc. Fourth Conf. Num. Sim. Plasmas, Vol. 3, 1970, p. 67. 



29 



